#LOAD packages
library("dummies")
## dummies-1.5.6 provided by Decision Patterns
library("AER")
## Loading required package: car
## Loading required package: carData
## Loading required package: lmtest
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: survival
library("plotly")
## Loading required package: ggplot2
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library('RColorBrewer')
library("data.table")
library("mlogit")
## Loading required package: Formula
library("gmnl")
## Loading required package: maxLik
## Loading required package: miscTools
##
## Please cite the 'maxLik' package as:
## Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
##
## If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
## https://r-forge.r-project.org/projects/maxlik/
####Part One --- Logit Model without segmentation#####
#Estimate a multinomial logit model using gmnl and mlogit
#Read data
rm(list=ls())
#getwd()
data=fread("2020 Spring A/MKT440 Pricing Analytics/project 2/kiwi_bubbles_P2.csv",stringsAsFactors = F)
#Load demographic data
demo=fread("2020 Spring A/MKT440 Pricing Analytics/project 2/demo_P2.csv",stringsAsFactors = F)
#Number of individuals
#Number of individuals
n = 1000
#Drop observations with stockout
data=data[!(data$price.KB==99),]
data=data[!(data$price.KR==99),]
data=data[!(data$price.MB==99),]
#Calculate elasticities using average price
avgPriceKB = mean(data$price.KB) #1.381332
avgPriceKR = mean(data$price.KR) #1.378087
avgPriceMB = mean(data$price.MB) #1.345585
avgPriceAll=c(avgPriceKB,avgPriceKR,avgPriceMB)
##################################### Functions we will use later#####################################
#Estimate multinomial logit model
#Convert data using 'mlogit.data'
mlogitdata = mlogit.data(data, id = 'id', varying = 4:7,choice = 'choice', shape = 'wide')
#Run MLE
mle = gmnl(choice ~ price, data = mlogitdata)
summary(mle)
##
## Model estimated on: Thu Feb 13 00:10:40 2020
##
## Call:
## gmnl(formula = choice ~ price, data = mlogitdata, method = "nr")
##
## Frequencies of categories:
##
## 0 KB KR MB
## 0.41564 0.18035 0.20039 0.20362
##
## The estimation took: 0h:0m:0s
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## KB:(intercept) 4.25316 0.32821 12.959 < 2.2e-16 ***
## KR:(intercept) 4.36240 0.32945 13.241 < 2.2e-16 ***
## MB:(intercept) 4.20440 0.31331 13.419 < 2.2e-16 ***
## price -3.73793 0.23671 -15.791 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Optimization of log-likelihood by Newton-Raphson maximisation
## Log Likelihood: -1909
## Number of observations: 1547
## Number of iterations: 4
## Exit of MLE: gradient close to zero
#beta0KB = 4.25316, beta0KR = 4.36240, beta0MB = 4.20440, beta1 = -3.73793
#beta0 of three products are similar, meaning they are selected equally frequently
beta1=-3.73793
beta0KB = 4.25316
beta0KR = 4.36240
beta0MB = 4.20440
#Functions:
#Demand:
#two proudcts:
para=c(beta0KB,beta0KR,beta0MB,beta1)
demand2=function(price1,price2,para){ #2 products, para contain all 4 numbers
prob=exp(para[1]+para[4]*price1)/(1+exp(para[1]+para[4]*price1)+exp(para[2]+para[4]*price2))
return(prob)
}
#three products:
demand1=function(priceKB,priceKR,priceMB,para){
prob=exp(para[1]+para[4]*priceKB)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
return(prob)
}
#agg_choice for two products
agg_choice2=function(priceKB,priceKR, own, rival) { #2 products
m=c(own+1,rival+1,5)
agg_choice=seg.share[1]*demand2(priceKB,priceKR,as.numeric(coef.est[1,m]))+
seg.share[2]*demand2(priceKB,priceKR,as.numeric(coef.est[2,m]))+
seg.share[3]*demand2(priceKB,priceKR,as.numeric(coef.est[3,m]))+
seg.share[4]*demand2(priceKB,priceKR,as.numeric(coef.est[4,m]))+
seg.share[5]*demand2(priceKB,priceKR,as.numeric(coef.est[5,m]))+
seg.share[6]*demand2(priceKB,priceKR,as.numeric(coef.est[6,m]))+
seg.share[7]*demand2(priceKB,priceKR,as.numeric(coef.est[7,m]))+
seg.share[8]*demand2(priceKB,priceKR,as.numeric(coef.est[8,m]))+
seg.share[9]*demand2(priceKB,priceKR,as.numeric(coef.est[9,m]))
return(agg_choice)
}
#agg_choice for three products
agg_choice=function(priceKB,priceKR,priceMB) {
agg_choice=seg.share[1]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[1,2:5]))+
seg.share[2]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[2,2:5]))+
seg.share[3]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[3,2:5]))+
seg.share[4]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[4,2:5]))+
seg.share[5]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[5,2:5]))+
seg.share[6]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[6,2:5]))+
seg.share[7]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[7,2:5]))+
seg.share[8]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[8,2:5]))+
seg.share[9]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[9,2:5]))
return(agg_choice)
}
#sum demand of KB+KR (based on 3 products environment)
demand_KB_KR=function(priceKB,priceKR,priceMB,para){
probKB=exp(para[1]+para[4]*priceKB)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
probKR=exp(para[2]+para[4]*priceKR)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
return(cbind(probKB,probKR))
}
demand_KB_KR_MB=function(priceKB,priceKR,priceMB,para){
probKB=exp(para[1]+para[4]*priceKB)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
probKR=exp(para[2]+para[4]*priceKR)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
probMB=exp(para[3]+para[4]*priceKR)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
return(cbind(probKB,probKR,probMB))
}
#sum profit of KB+KR (based on 3 products environment)
profit_KB_KR=function(priceKB,priceKR, priceMB,para){
profitKB=1000*(demand_KB_KR(priceKB,priceKR, priceMB,para)[,1]*(priceKB-0.5))
profitKR=1000*(demand_KB_KR(priceKB,priceKR, priceMB,para)[,2]*(priceKR-0.5))
return(cbind(profitKB,profitKR))
}
#################Part 3##############
###Elasticity
#beta1=-3.73793
#beta0KB = 4.25316
#beta0KR = 4.36240
#beta0MB = 4.20440
avgPriceKB = mean(data$price.KB) #1.381332
avgPriceKR = mean(data$price.KR) #1.378087
avgPriceMB = mean(data$price.MB) #1.345585
#KB+KR+MB
demandKB = demand1(avgPriceKB,avgPriceKR,avgPriceMB,para)
demandKR = demand1(avgPriceKR,avgPriceKB,avgPriceMB,para[c(2,1,3,4)])
demandMB = demand1(avgPriceMB,avgPriceKB,avgPriceKR,para[c(3,1,2,4)])
#Own-price elasticity of KB
ownElasKB = -(beta1)*avgPriceKB*(1-demandKB) #4.222673
#Own-price elasticity of KR
ownElasKR = -(beta1)*avgPriceKR*(1-demandKR) #4.143605
#Own-price elasticity of MB
ownElasMB = -(beta1)*avgPriceMB*(1-demandMB) #4.072252
#Cross-elasticities
crossElasKR = -(beta1)*avgPriceKR*demandKR #1.012787 -> canibalization
crossElasMB = -(beta1)*avgPriceMB*demandMB #0.9574505
crossElasKB = -(beta1)*avgPriceKB*demandKB #0.9188917
#CREATE DATA FRAME
elas_KB_KR_MB <- data.frame("Product"=c("KB","KR","MB"), "Own Elasticity"=c(ownElasKB,ownElasKR,ownElasMB),
"Cross Price Elasticity"=c(crossElasKB,crossElasKR,crossElasMB))
elas_KB_KR_MB
## Product Own.Elasticity Cross.Price.Elasticity
## 1 KB 4.257845 0.9054761
## 2 KR 4.131272 1.0199190
## 3 MB 4.069542 0.9601601
############2. Optimal Price for KB AND KR
#Multinomial logit: illustration
#Unit cost
uc=0.5
#KB:(intercept) KR:(intercept) MB:(intercept) price
#4.253157 4.362403 4.204396 -3.737931
#Set parameter
#The first element of "para" is beta0KB, the second element is beta0KR and the third beta1, FORTH IS beta0MB.
priceMB=1.43
demand_KB_KR=function(priceKB,priceKR,priceMB,para){
probKB=exp(para[1]+para[4]*priceKB)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
probKR=exp(para[2]+para[4]*priceKR)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
return(cbind(probKB,probKR))
}
#Write profit as a function of prices we set and model parameters
profit_KB_KR=function(priceKB,priceKR, priceMB,para){
profitKB=1000*(demand_KB_KR(priceKB,priceKR, priceMB,para)[,1]*(priceKB-0.5))
profitKR=1000*(demand_KB_KR(priceKB,priceKR, priceMB,para)[,2]*(priceKR-0.5))
return(cbind(profitKB,profitKR))
}
#Choose space of prices to search for the optimal price over
#Because we search over two dimensions, create complete combination of the two prices
#Choose space of prices to search for the optimal price over
aux=seq(1,3,0.01)
#Because we search over two dimensions, create complete combination
#of the two prices
pricespace=expand.grid(aux,aux)
#Compute profit at each realization of this price space.
#I write for-loop. While there are ways to get this done in a vectorized way,
#this for-loop code probably helps some in project 2.
#Calculate profit
#At each iteration of the loop, I take one realization of [P^KB,P^KR] pair and evaluate
#profit at that realization.
profitmat=matrix(0L,nrow(pricespace),1)
for (i in 1:nrow(pricespace)){
profitmat[i]=sum(profit_KB_KR(pricespace[i,1],pricespace[i,2],1.43,para))
}
#Draw figure
xaxis=list(title="P^{KB}")
yaxis=list(autorange = "reversed",title="P^{KR}")
zaxis=list(title="Profit")
p=plot_ly(x=pricespace[,1],y=pricespace[,2],z=as.numeric(profitmat),
type="scatter3d",mode="markers",
marker = list(color = as.numeric(profitmat), colorscale = c('#FFE1A1', '#683531'), showscale = TRUE))%>%
layout(scene=list(xaxis=xaxis,yaxis=yaxis,zaxis=zaxis))%>%
config(mathjax = 'cdn')
p
max(profitmat) #393.408; optimal price for KB is 1.16, optimal price for KR is 1.16
## [1] 393.408
pricespace[profitmat==max(profitmat)] #optimal price for KB is 1.16, optimal price for KR is 1.16
## [1] 1.16 1.16
#######Part 4######################
set.seed(0)
###### try and find the number of segments
#Clustering
#number of individuals in demo
N=283
#Try optimal Number of Clusters
#Clustering&kmeans
#cluster=8
segList = list()
#create a copy of original data
data2 = data
#write a for-loop to test what is the reasonable number of segments
for(i in 1:18){
demo_cluster = kmeans(x=demo[, 2:18], centers = i, nstart = 1000)
data2=data
# now combine cluster identity into the raw data
cluster_id = data.frame(id = demo$id)
cluster_id$cluster = demo_cluster$cluster
data2 = merge(data2, cluster_id, by = "id", all.x = T)
# for those who don't fit in any cluster, group them into one additional cluster
data2$cluster[is.na(data2$cluster)] = i+1
#segment share
seg.share = c( table(demo_cluster$cluster),N - sum(table(demo_cluster$cluster))) / N
#input seg.share into segTable
segList[[i]] = seg.share
#recover original data
data2 = data2[,1:8]
}
segList[[8]]
## 1 2 3 4 5 6 7
## 0.10954064 0.16607774 0.10954064 0.14134276 0.12014134 0.13780919 0.13427562
## 8
## 0.08127208 0.00000000
segList[[9]]
## 1 2 3 4 5 6 7
## 0.10954064 0.01413428 0.11307420 0.13427562 0.14840989 0.07067138 0.16961131
## 8 9
## 0.13780919 0.10247350 0.00000000
#we can see from the segList that when we set the number of segments to 9, the segment share becomes very small (0.014)
#which means that segment contains less than 4 people
#Therefore, we choose 8 segments
##############use 8 segments
demo_cluster = kmeans(x=demo[, 2:18], centers = 8, nstart = 1000)
summary(demo_cluster)
## Length Class Mode
## cluster 283 -none- numeric
## centers 136 -none- numeric
## totss 1 -none- numeric
## withinss 8 -none- numeric
## tot.withinss 1 -none- numeric
## betweenss 1 -none- numeric
## size 8 -none- numeric
## iter 1 -none- numeric
## ifault 1 -none- numeric
# now combine cluster identity into the raw data
cluster_id = data.frame(id = demo$id)
cluster_id$cluster = demo_cluster$cluster
data = merge(data, cluster_id, by = "id", all.x = T)
# for those who don't fit in any cluster, group them into one additional cluster
data$cluster[is.na(data$cluster)] = 9
#segment share
seg.share = c( table(demo_cluster$cluster),N - sum(table(demo_cluster$cluster))) / N
#Estimate multinomial logit model separately for each segment
#store the coefficients
coef.est = data.frame(segment = 1:9, intercept.KB = NA, intercept.KR = NA,
intercept.MB = NA, price.coef = NA)
for (seg in 1:9) {
# During each loop, pick subset of data of consumers from each segment.
data.sub = subset(data, cluster == seg)
#Using that data, the rest remains the same.
mlogitdata = mlogit.data(data.sub,id="id",varying=4:7,choice="choice",shape="wide")
#Run MLE.
mle = gmnl(choice ~ price, data = mlogitdata)
mle
#Store the outcome in the coef.est matrix.
coef.est[seg, 2:5] = mle$coefficients
}
coef.est
## segment intercept.KB intercept.KR intercept.MB price.coef
## 1 1 3.8689983 4.354056 4.0523865 -3.502896
## 2 2 3.9972969 3.958938 3.8837551 -3.715366
## 3 3 2.9694016 3.867674 2.7276100 -2.909001
## 4 4 7.3034174 7.138563 7.1181389 -5.793619
## 5 5 2.3336806 3.112574 2.9252012 -2.896447
## 6 6 7.6063946 6.661934 7.4775392 -5.897474
## 7 7 4.8086045 4.606528 5.6294806 -4.517302
## 8 8 0.9169255 1.673183 0.4573439 -1.251711
## 9 9 5.1174301 4.509340 4.5449628 -4.062526
popularityCompare = data.frame(segments = 1:9, differKR = coef.est$intercept.KR-coef.est$intercept.KB,
differMB = coef.est$intercept.MB-coef.est$intercept.KB)
popularityCompare
## segments differKR differMB
## 1 1 0.48505759 0.1833882
## 2 2 -0.03835883 -0.1135418
## 3 3 0.89827253 -0.2417916
## 4 4 -0.16485417 -0.1852785
## 5 5 0.77889343 0.5915207
## 6 6 -0.94446055 -0.1288554
## 7 7 -0.20207671 0.8208761
## 8 8 0.75625710 -0.4595815
## 9 9 -0.60808976 -0.5724673
differCompare = data.frame(differKR = mean(popularityCompare$differKR), differMB = mean(popularityCompare$differMB))
differCompare
## differKR differMB
## 1 0.1067378 -0.0117479
################################################################
######Elasticity
############### ############### PART (4)################### ###############
############### ############### PART (4)################### ###############
#AVG PRICE (3 products, kb kr mb) for all 3 segments
demand1=function(priceKB,priceKR,priceMB,para){
probKB=exp(para[1]+para[4]*priceKB)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
probKR=exp(para[2]+para[4]*priceKR)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
probMB=exp(para[3]+para[4]*priceKR)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
return(cbind(probKB,probKR,probMB))
}
agg_choice=function(priceKB,priceKR,priceMB) { #3 products segment
agg_choice=seg.share[1]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[1,2:5]))+
seg.share[2]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[2,2:5]))+
seg.share[3]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[3,2:5]))+
seg.share[4]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[4,2:5]))+
seg.share[5]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[5,2:5]))+
seg.share[6]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[6,2:5]))+
seg.share[7]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[7,2:5]))+
seg.share[8]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[8,2:5]))+
seg.share[9]*demand1(priceKB,priceKR,priceMB,as.numeric(coef.est[9,2:5]))
return(agg_choice)
}
avgPriceVec=list()
for (i in 1:3){
avgPriceKB = mean(data$price.KB[data$cluster==i])
avgPriceKR = mean(data$price.KR[data$cluster==i])
avgPriceMB = mean(data$price.MB[data$cluster==i])
avgPriceVec[[i]]=c(avgPriceKB,avgPriceKR,avgPriceMB)
}
###aggregate_choice and demand#########
para2=c(beta0KR,beta0MB,beta1)
demand2=function(price1,price2,para2){ #2 products, para contains 3 #s, func SAME NAME
prob1=exp(para2[1]+para2[3]*price1)/(1+exp(para2[1]+para2[3]*price1)+exp(para2[2]+para2[3]*price2))
prob2=exp(para2[2]+para2[3]*price1)/(1+exp(para2[1]+para2[3]*price1)+exp(para2[2]+para2[3]*price2))
return(cbind(prob1,prob2))
}
agg_choice2=function(priceKB,priceKR, own, rival) { #2 products segment
m=c(own+1,rival+1,5)
agg_choice=seg.share[1]*demand2(priceKB,priceKR,as.numeric(coef.est[1,m]))+
seg.share[2]*demand2(priceKB,priceKR,as.numeric(coef.est[2,m]))+
seg.share[3]*demand2(priceKB,priceKR,as.numeric(coef.est[3,m]))+
seg.share[4]*demand2(priceKB,priceKR,as.numeric(coef.est[4,m]))+
seg.share[5]*demand2(priceKB,priceKR,as.numeric(coef.est[5,m]))+
seg.share[6]*demand2(priceKB,priceKR,as.numeric(coef.est[6,m]))+
seg.share[7]*demand2(priceKB,priceKR,as.numeric(coef.est[7,m]))+
seg.share[8]*demand2(priceKB,priceKR,as.numeric(coef.est[8,m]))+
seg.share[9]*demand2(priceKB,priceKR,as.numeric(coef.est[9,m]))
return(agg_choice)
}
#own elasticity
elasticity_own=function(avgPriceAll,own){
coef2=c(1,2,3)[c(1,2,3)!=own][1]
coef3=c(1,2,3)[c(1,2,3)!=own][2]
agg_choice=agg_choice(avgPriceAll[own],avgPriceAll[coef2],avgPriceAll[coef3])
agg_choicenew=agg_choice(avgPriceAll[own]*1.01,avgPriceAll[coef2],avgPriceAll[coef3])
(agg_choicenew-agg_choice)/agg_choice
}
KB4=(elasticity_own(avgPriceAll,1)[1])*-100 #KB -0.04392695
KR4=(elasticity_own(avgPriceAll,2)[1])*-100 #KR -0.04356302
MB4=(elasticity_own(avgPriceAll,3)[1])*-100 #MB -0.04112575
Own_elas_KB_KR_MB <- data.frame("Product"=c("KB","KR","MB"), "Own Elasticity"=c(KB4,KR4,MB4))
Own_elas_KB_KR_MB
## Product Own.Elasticity
## 1 KB 4.393557
## 2 KR 4.375514
## 3 MB 4.186238
#cross elasticity
elasticity_cross=function(avgPriceAll,own){
e=vector()
coef2=c(1,2,3)[c(1,2,3)!=own][1]
coef3=c(1,2,3)[c(1,2,3)!=own][2]
agg_choice=agg_choice(avgPriceAll[own],avgPriceAll[coef2],avgPriceAll[coef3])[1]
agg_choicenew1=agg_choice(avgPriceAll[own],avgPriceAll[coef2]*1.01,avgPriceAll[coef3])[1]
agg_choicenew2=agg_choice(avgPriceAll[own],avgPriceAll[coef2],avgPriceAll[coef3]*1.01)[1]
e=append(e, (agg_choicenew1-agg_choice)/agg_choice)
e=append(e, (agg_choicenew2-agg_choice)/agg_choice)
e
}
Cros_KB_KR_MB=(elasticity_cross(avgPriceAll, 1))*100
Cros_KB_KR_MB
## [1] 0.9459524 1.1224485
Cros_KR_KB_MB=(elasticity_cross(avgPriceAll, 2))*100
Cros_KR_KB_MB
## [1] 0.9363006 1.1235926
Cros_MB_KB_KR=(elasticity_cross(avgPriceAll, 3))*100
Cros_MB_KB_KR
## [1] 0.9351753 1.0089358
###### (4) last question###########
#Profit
#given MB, look at KR
pricespace=seq(1,3,0.01)
uc=0.5
profit=1000*agg_choice2(pricespace,1.43,2,3)[,1]*(pricespace-uc)
priceKRBest=pricespace[profit==max(profit)]
priceKRBest#best price for KR=1.07
## [1] 1.07
max(profit) #max profit for KR 294.4419
## [1] 294.4419
#profit for MB
profit_MB=1000*agg_choice2(1.43,priceKRBest,3,2)[,1]*(1.43-uc) #109.1702
profit_MB
## prob1
## 109.1702
#GIVEN MB, look at KR, KB
para=c(beta0KB,beta0KR,beta0MB,beta1)
profit_KB_KR_seg=function(priceKB,priceKR, priceMB,para){
profitKB=1000*(agg_choice(priceKB,priceKR, priceMB)[,1]*(priceKB-0.5))
profitKR=1000*(agg_choice(priceKB,priceKR, priceMB)[,2]*(priceKR-0.5))
return(cbind(profitKB,profitKR))
}
aux=seq(1,3,0.01)
#Because we search over two dimensions, create complete combination
#of the two prices
pricespace=expand.grid(aux,aux)
profitmat=matrix(0L,nrow(pricespace),1)
for (i in 1:nrow(pricespace)){
profitmat[i]=sum(profit_KB_KR_seg(pricespace[i,1],pricespace[i,2],1.43,para))
}
max(profitmat) #387.6119;
## [1] 387.6119
priceKB_KRBest=pricespace[profitmat==max(profitmat)] #optimal price for KB is 1.15, optimal price for KR is 1.19
priceKB_KRBest
## [1] 1.15 1.19
#MB new profit
profit_MB2=1000*agg_choice(1.43,priceKB_KRBest,priceKB_KRBest)[2]*(1.43-uc) #244.3255 new MB profit
profit_MB2 #89.63174
## [1] 89.63174
#part 5
# round 1
#MANGO:
pricespace=seq(0,2,0.01)
uc=0.5
profit=1000*agg_choice(pricespace,1.15,1.19)[,1]*(pricespace-uc)
max(profit) #170.9814
## [1] 170.9814
pricespace[profit==max(profit)] #0.96
## [1] 0.96
#KB KR
#Choose space of prices to search for the optimal price over
aux=seq(1,3,0.01)
#Because we search over two dimensions, create complete combination
#of the two prices
pricespace=expand.grid(aux,aux)
profitmat=matrix(0L,nrow(pricespace),1)
for (i in 1:nrow(pricespace)){
profitmat[i]=sum(profit_KB_KR_seg(pricespace[i,1],pricespace[i,2],0.96,para))
}
max(profitmat) #269.6908
## [1] 269.6908
pricespace[profitmat==max(profitmat)] #1.01 1.09
## [1] 1.01 1.09
####Round 2
#MANGO:
pricespace=seq(0,2,0.01)
uc=0.5
profit=1000*agg_choice(pricespace,1.01,1.09)[,1]*(pricespace-uc)
max(profit) #140.0129
## [1] 140.0129
pricespace[profit==max(profit)] #0.92
## [1] 0.92
#KB KR
#Choose space of prices to search for the optimal price over
aux=seq(1,3,0.01)
#Because we search over two dimensions, create complete combination
#of the two prices
pricespace=expand.grid(aux,aux)
profitmat=matrix(0L,nrow(pricespace),1)
for (i in 1:nrow(pricespace)){
profitmat[i]=sum(profit_KB_KR_seg(pricespace[i,1],pricespace[i,2],0.92,para))
}
max(profitmat) #256.2545
## [1] 256.2545
pricespace[profitmat==max(profitmat)] #1.00 1.08
## [1] 1.00 1.08
####Round 3
#MANGO:
pricespace=seq(0,2,0.01)
uc=0.5
profit=1000*agg_choice(pricespace,1.00,1.08)[,1]*(pricespace-uc)
max(profit) #137.4193
## [1] 137.4193
pricespace[profit==max(profit)] #0.92
## [1] 0.92
#KB KR
#Choose space of prices to search for the optimal price over
aux=seq(1,3,0.01)
#Because we search over two dimensions, create complete combination
#of the two prices
pricespace=expand.grid(aux,aux)
profitmat=matrix(0L,nrow(pricespace),1)
for (i in 1:nrow(pricespace)){
profitmat[i]=sum(profit_KB_KR_seg(pricespace[i,1],pricespace[i,2],0.92,para))
}
max(profitmat) #256.2545
## [1] 256.2545
pricespace[profitmat==max(profitmat)] #1.00 1.08
## [1] 1.00 1.08